#R version 4.2.1 (2022-06-23)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 10 x64 (build 19044)

library(dplyr) #version 1.0.9
library(readr) #version 2.1.2
library(stringr) #version 1.4.0
library(zoo) #version 1.8-10
library(sandwich) #version 3.0-2
library(lmtest) #version 0.9-40
library(prediction) #version 0.3.14
library(ggpubr) #version 0.4.0
library(estimatr) #version 1.0.0
library(margins) #version 0.3.26
library(broom) #version 1.0.0
library(scales) #version 1.2.0
library(stargazer) #version 5.2.3

`%notin%` <- function(x,y) !(x %in% y) 

#Table 6 - GWF Additive and Multiplicative Models####
#institutionalization
d1 <- read_csv("data/gwf_fa_pi.csv") 
d2 <- read_csv("data/gwf_irt.csv") %>% 
  rename(gwf_casename = case) %>%
  select(gwf_casename, transition, gwf_irt) 

gwf1 <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("gwf_casename", "year")) %>% 
  left_join(d2, by = c("gwf_casename")) %>% 
  group_by(gwf_casename) %>% 
  mutate_at(.vars = vars(gwf_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, gwf_casename, gwf_regimetype, v2xps_party, v2pssunpar, gwf_incumbent_pi, gwf_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, gwf_priorduration, coups, gwf_number, EFindex, theta_mean) %>%
  mutate(gwf_incumbent_pi = ifelse(is.na(gwf_incumbent_pi), v2xps_party, gwf_incumbent_pi)) %>% 
  rename(incumbent_pi = gwf_incumbent_pi) %>% 
  filter(transition == year, gwf_casename %notin% c("Bolivia 43-46", "Thailand 57-73", "Pakistan 58-71", "Haiti 88-90",
                                                    "Sudan 58-64", "Ecuador 63-66", "Venezuela 48-58")) %>% 
  ungroup() %>% 
  unique()

#strength
d1 <- read_csv("data/gwf_fa_strength.csv") 
d2 <- read_csv("data/gwf_irt.csv") %>% 
  rename(gwf_casename = case) %>%
  dplyr::select(gwf_casename, transition, gwf_irt) 

gwf2 <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("gwf_casename", "year")) %>% 
  left_join(d2, by = c("gwf_casename")) %>% 
  group_by(gwf_casename) %>% 
  mutate_at(.vars = vars(gwf_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

gwf2$v2pssunpar <- scales::rescale(gwf2$v2pssunpar, to = c(0, 1), from = range(gwf2$v2pssunpar, na.rm = TRUE, finite = TRUE))

gwf2 %<>%
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, gwf_casename, gwf_regimetype, v2xps_party, v2pssunpar, gwf_incumbent_strength, 
                gwf_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, gwf_priorduration, coups, gwf_number, EFindex, theta_mean) %>%
  mutate(gwf_incumbent_strength = ifelse(is.na(gwf_incumbent_strength), v2pssunpar, gwf_incumbent_strength)) %>% 
  rename(incumbent_strength = gwf_incumbent_strength) %>% 
  filter(transition == year, gwf_casename %notin% c("Bolivia 43-46", "Thailand 57-73", "Pakistan 58-71", "Haiti 88-90",
                                                    "Sudan 58-64", "Ecuador 63-66", "Venezuela 48-58")) %>% 
  ungroup() %>% 
  unique() %>% 
  select(gwf_casename, incumbent_strength)

gwf <- left_join(gwf1, gwf2, by = "gwf_casename")

m1 <- lm(gwf_irt ~ incumbent_pi + incumbent_strength, data = gwf) 

m2 <- lm(gwf_irt ~ incumbent_pi + incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + gwf_priorduration + coups + gwf_number + EFindex + theta_mean, data = gwf) 

m3 <- lm(gwf_irt ~ incumbent_pi + incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + gwf_priorduration + coups + gwf_number + EFindex + theta_mean, data = gwf) 

m4 <- lm(gwf_irt ~ incumbent_pi*incumbent_strength, data = gwf) 

m5 <- lm(gwf_irt ~ incumbent_pi*incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + gwf_priorduration + coups + gwf_number + EFindex + theta_mean, data = gwf) 

m6 <- lm(gwf_irt ~ incumbent_pi*incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + gwf_priorduration + coups + gwf_number + EFindex + theta_mean, data = gwf) 

stargazer(m1, m2, m3, m4, m5, m6, se = starprep(m1, m2, m3, m4, m5, m6, se_type = "HC1"), omit.stat=c("f", "ser"), omit = "Constant",
          order = c(1, 2, 18),  covariate.labels = c("Party Institutionalization", "Party Strength", "Institutionalization*Strength", "GINI", "Rural Inequality", 
                                                     "Military Spending", "Area", "Urban Pop", "GDP PC", 
                                                     "Resource Wealth", "Islam", "Region", "Education", "Duration", "Coups", "Regimes", "Ethnic Frac", "Repression"),
          column.labels   = c("Additive", "Multiplicative"),
          column.separate = c(3, 3))



#Table 7 - PAR Additive and Multiplicative Models####
#institutionalization
d1 <- read_csv("data/svolik_fa_pi.csv") 
d2 <- read_csv("data/svolik_irt.csv") %>% 
  rename(svolik_case = case) %>%
  select(svolik_case, transition, svolik_irt)

svolik1 <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("svolik_case", "year")) %>% 
  left_join(d2, by = c("svolik_case")) %>% 
  group_by(svolik_case) %>% 
  mutate_at(.vars = vars(svolik_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, svolik_case, svolik_regime, v2xps_party, v2pssunpar, svolik_incumbent_pi, svolik_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea,
         svolik_priorduration, coups, svolik_number, EFindex, theta_mean) %>%
  mutate(svolik_incumbent_pi = ifelse(is.na(svolik_incumbent_pi), v2xps_party, svolik_incumbent_pi)) %>% 
  rename(incumbent_pi = svolik_incumbent_pi)  %>% 
  filter(transition == year, svolik_case %notin% c("Bolivia 1946-1946", "Thailand 1971-1972", "Pakistan 1958-1972", "Haiti 1986-1990", 
                                                   "Sudan 1959-1964", "Ecuador 1963-1965", "Haiti 1950-1956",
                                                   "Comoros 1999-2001", "Venezuela 1948-1949", "Burkina Faso (Upper Volta) 1966-1971", "Peru 1949-1956"))%>% 
  ungroup() %>% 
  unique()

#strength
d1 <- read_csv("data/svolik_fa_strength.csv") 
d2 <- read_csv("data/svolik_irt.csv") %>% 
  rename(svolik_case = case) %>%
  dplyr::select(svolik_case, transition, svolik_irt)

svolik2 <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("svolik_case", "year")) %>% 
  left_join(d2, by = c("svolik_case")) %>% 
  group_by(svolik_case) %>% 
  mutate_at(.vars = vars(svolik_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

svolik2$v2pssunpar <- scales::rescale(svolik2$v2pssunpar, to = c(0, 1), from = range(svolik2$v2pssunpar, na.rm = TRUE, finite = TRUE))

svolik2 %<>%  
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar))%>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, svolik_case, svolik_regime, v2xps_party, v2pssunpar, svolik_incumbent_strength, 
                svolik_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea,
                svolik_priorduration, coups, svolik_number, EFindex, theta_mean) %>%
  mutate(svolik_incumbent_strength = ifelse(is.na(svolik_incumbent_strength), v2pssunpar, svolik_incumbent_strength)) %>% 
  rename(incumbent_strength = svolik_incumbent_strength) %>% 
  filter(transition == year, svolik_case %notin% c("Bolivia 1946-1946", "Thailand 1971-1972", "Pakistan 1958-1972", "Haiti 1986-1990", 
                                                   "Sudan 1959-1964", "Ecuador 1963-1965", "Haiti 1950-1956",
                                                   "Comoros 1999-2001", "Venezuela 1948-1949", "Burkina Faso (Upper Volta) 1966-1971", "Peru 1949-1956"))%>% 
  ungroup() %>% 
  unique()%>% 
  select(svolik_case, incumbent_strength)

svolik <- left_join(svolik1, svolik2, by = "svolik_case")


m1 <- lm(svolik_irt ~ incumbent_pi + incumbent_strength , data = svolik) 

m2 <- lm(svolik_irt ~ incumbent_pi + incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik) 

m3 <- lm(svolik_irt ~ incumbent_pi + incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik) 

m4 <- lm(svolik_irt ~ incumbent_pi*incumbent_strength, data = svolik)

m5 <- lm(svolik_irt ~ incumbent_pi*incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik)

m6 <- lm(svolik_irt ~ incumbent_pi*incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik) 

stargazer(m1, m2, m3, m4, m5, m6, se = starprep(m1, m2, m3, m4, m5, m6, se_type = "HC1"), omit.stat=c("f", "ser"), omit = "Constant",
          order = c(1, 2, 18),  covariate.labels = c("Party Institutionalization", "Party Strength", "Institutionalization*Strength", "GINI", "Rural Inequality", 
                                                     "Military Spending", "Area", "Urban Pop", "GDP PC", 
                                                     "Resource Wealth", "Islam", "Region", "Education", "Duration", "Coups", "Regimes", "Ethnic Frac", "Repression"),
          column.labels   = c("Additive", "Multiplicative"),
          column.separate = c(3, 3))




#Table 8 - DD Additive and Multiplicative Models#####
#institutionalization
d1 <- read_csv("data/dd_fa_pi.csv") 
d2 <- read_csv("data/dd_irt.csv") %>% 
  rename(dd_case = case) %>%
  select(dd_case, transition, dd_irt)

dd1 <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("dd_case", "year")) %>% 
  left_join(d2, by = c("dd_case")) %>% 
  group_by(dd_case) %>% 
  mutate_at(.vars = vars(dd_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, dd_case, dd_regime, v2xps_party, v2pssunpar, dd_incumbent_pi, dd_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
         dd_priorduration, coups, dd_number, EFindex, theta_mean) %>%
  mutate(dd_incumbent_pi = ifelse(is.na(dd_incumbent_pi), v2xps_party, dd_incumbent_pi)) %>% 
  rename(incumbent_pi = dd_incumbent_pi)  %>% 
  filter(transition == year, dd_case %notin% c("Bolivia 1946-1946", "Thailand 1946-1972", "Pakistan 1958-1970", "Haiti 1986-1989", 
                                               "Sudan 1958-1963", "Ecuador 1963-1965", "Haiti 1950-1955",
                                               "Comoros 1999-2003", "Venezuela 1948-1958")) %>% 
  ungroup() %>% 
  unique()

#strength
d1 <- read_csv("data/dd_fa_strength.csv") 
d2 <- read_csv("data/dd_irt.csv") %>% 
  rename(dd_case = case) %>%
  dplyr::select(dd_case, transition, dd_irt)

dd2 <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("dd_case", "year")) %>% 
  left_join(d2, by = c("dd_case")) %>% 
  group_by(dd_case) %>% 
  mutate_at(.vars = vars(dd_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

dd2$v2pssunpar <- scales::rescale(dd2$v2pssunpar, to = c(0, 1), from = range(dd2$v2pssunpar, na.rm = TRUE, finite = TRUE))

dd2 %<>%
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, dd_case, dd_regime, v2xps_party, v2pssunpar, dd_incumbent_strength, 
                dd_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
                dd_priorduration, coups, dd_number, EFindex, theta_mean) %>%
  mutate(dd_incumbent_strength = ifelse(is.na(dd_incumbent_strength), v2pssunpar, dd_incumbent_strength)) %>% 
  rename(incumbent_strength = dd_incumbent_strength) %>% 
  filter(transition == year, dd_case %notin% c("Bolivia 1946-1946", "Thailand 1946-1972", "Pakistan 1958-1970", "Haiti 1986-1989", 
                                               "Sudan 1958-1963", "Ecuador 1963-1965", "Haiti 1950-1955",
                                               "Comoros 1999-2003", "Venezuela 1948-1958")) %>% 
  ungroup() %>% 
  unique()  %>% 
  select(dd_case, incumbent_strength)

dd <- left_join(dd1, dd2, by = "dd_case")

m1 <- lm(dd_irt ~ incumbent_pi + incumbent_strength, data = dd) 

m2 <- lm(dd_irt ~ incumbent_pi + incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + dd_priorduration + coups + dd_number + EFindex + theta_mean, data = dd) 

m3 <- lm(dd_irt ~ incumbent_pi + incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + dd_priorduration + coups + dd_number + EFindex + theta_mean, data = dd) 

m4 <- lm(dd_irt ~ incumbent_pi*incumbent_strength, data = dd)

m5 <- lm(dd_irt ~ incumbent_pi*incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + dd_priorduration + coups + dd_number + EFindex + theta_mean, data = dd)

m6 <- lm(dd_irt ~ incumbent_pi*incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + dd_priorduration + coups + dd_number + EFindex + theta_mean, data = dd) 

stargazer(m1, m2, m3, m4, m5, m6, se = starprep(m1, m2, m3, m4, m5, m6, se_type = "HC1"), omit.stat=c("f", "ser"), omit = "Constant",
          order = c(1, 2, 18),  covariate.labels = c("Party Institutionalization", "Party Strength", "Institutionalization*Strength", "GINI", "Rural Inequality", 
                                                     "Military Spending", "Area", "Urban Pop", "GDP PC", 
                                                     "Resource Wealth", "Islam", "Region", "Education", "Duration", "Coups", "Regimes", "Ethnic Frac", "Repression"),
          column.labels   = c("Additive", "Multiplicative"),
          column.separate = c(3, 3))

#Table 9 - WTH Additive and Multiplicative Models#####
#institutionalization
d1 <- read_csv("data/wth_fa_pi.csv") 
d2 <- read_csv("data/wth_irt.csv") %>% 
  rename(wth_case = case) %>%
  select(wth_case, transition, wth_irt)

wth1 <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("wth_case", "year")) %>% 
  left_join(d2, by = c("wth_case")) %>% 
  group_by(wth_case) %>% 
  mutate_at(.vars = vars(wth_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, wth_case, wth_regime, v2xps_party, v2pssunpar, wth_incumbent_pi, wth_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
         wth_priorduration, coups, wth_number, EFindex, theta_mean) %>%
  mutate(wth_incumbent_pi = ifelse(is.na(wth_incumbent_pi), v2xps_party, wth_incumbent_pi)) %>% 
  rename(incumbent_pi = wth_incumbent_pi)  %>% 
  filter(transition == year, wth_case %notin% c("Haiti 1986-1989", "Comoros 1999-2001")) %>% 
  ungroup() %>% 
  unique()

#strength
d1 <- read_csv("data/wth_fa_strength.csv") 
d2 <- read_csv("data/wth_irt.csv") %>% 
  rename(wth_case = case) %>%
  dplyr::select(wth_case, transition, wth_irt)

wth2 <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("wth_case", "year")) %>% 
  left_join(d2, by = c("wth_case")) %>% 
  group_by(wth_case) %>% 
  mutate_at(.vars = vars(wth_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

wth2$v2pssunpar <- scales::rescale(wth2$v2pssunpar, to = c(0, 1), from = range(wth2$v2pssunpar, na.rm = TRUE, finite = TRUE))

wth2 %<>%
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, wth_case, wth_regime, v2xps_party, v2pssunpar, wth_incumbent_strength, 
                wth_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
                wth_priorduration, coups, wth_number, EFindex, theta_mean) %>%
  mutate(wth_incumbent_strength = ifelse(is.na(wth_incumbent_strength), v2pssunpar, wth_incumbent_strength)) %>% 
  rename(incumbent_strength = wth_incumbent_strength) %>% 
  filter(transition == year, wth_case %notin% c("Haiti 1986-1989", "Comoros 1999-2001")) %>% 
  ungroup() %>% 
  unique()  %>% 
  select(wth_case, incumbent_strength)

wth <- left_join(wth1, wth2, by = "wth_case")

m1 <- lm(wth_irt ~ incumbent_pi + incumbent_strength, data = wth) 

m2 <- lm(wth_irt ~ incumbent_pi + incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + wth_priorduration + coups + wth_number + EFindex + theta_mean, data = wth) 

m3 <- lm(wth_irt ~ incumbent_pi + incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + wth_priorduration + coups + wth_number + EFindex + theta_mean, data = wth) 

m4 <- lm(wth_irt ~ incumbent_pi*incumbent_strength, data = wth)

m5 <- lm(wth_irt ~ incumbent_pi*incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + wth_priorduration + coups + wth_number + EFindex + theta_mean, data = wth)

m6 <- lm(wth_irt ~ incumbent_pi*incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + wth_priorduration + coups + wth_number + EFindex + theta_mean, data = wth) 

stargazer(m1, m2, m3, m4, m5, m6, se = starprep(m1, m2, m3, m4, m5, m6, se_type = "HC1"), omit.stat=c("f", "ser"), omit = "Constant",
          order = c(1, 2, 18),  covariate.labels = c("Party Institutionalization", "Party Strength", "Institutionalization*Strength", "GINI", "Rural Inequality", 
                                                     "Military Spending", "Area", "Urban Pop", "GDP PC", 
                                                     "Resource Wealth", "Islam", "Region", "Education", "Duration", "Coups", "Regimes", "Ethnic Frac", "Repression"),
          column.labels   = c("Additive", "Multiplicative"),
          column.separate = c(3, 3))


#Table 5####

c1 <- cor(select(gwf, incumbent_pi, incumbent_strength))[1,2]
c2 <- cor(select(svolik, incumbent_pi, incumbent_strength))[1,2]
c3 <- cor(select(dd, incumbent_pi, incumbent_strength))[1,2]
c4 <- cor(select(wth, incumbent_pi, incumbent_strength))[1,2]

as.data.frame(c(c1, c2, c3, c4)) %>% 
  mutate(dataset = c("GWF", "PAR", "DD", "WTH"))

